Stability and Bifurcation Analysis of Coupled Fitzhugh-Nagumo Oscillators

Neurons are the central biological objects in understanding how the brain works. The famous Hodgkin-Huxley model, which describes how action potentials of a neuron are initiated and propagated, consists of four coupled nonlinear differential equations. Because these equations are difficult to deal with, there also exist several simplified models, of which many exhibit polynomial-like non-linearity. Examples of such models are the Fitzhugh-Nagumo (FHN) model, the Hindmarsh-Rose (HR) model, the MorrisLecar (ML) model and the Izhikevich model. In this work, we first prescribe the biologically relevant parameter ranges for the FHN model and subsequently study the dynamical behaviour of coupled neurons on small networks of two or three nodes. To do this, we use a computational real algebraic geometry method called the Discriminant Variety (DV) method to perform the stability and bifurcation analysis of these small networks. A time series analysis of the FHN model can be found elsewhere in related work [15]. The Fitzhugh-Nagumo Model The Fitzhugh-Nagumo (FHN) system of equations is a prototypic model of excitability. Introduced by FitzHugh [13] with equivalent circuit by Nagumo et al.[25] the system is a generalization of the Van der Pol oscillator [11] and a reduction of the neural electrophysiological model due to Hodgkin and Huxley (HH) [16]. As a generic model of excitability and oscillatory dynamical behavior, the FHN system is of relevance for a range of physical and physiological research topics[21, 27, 22, 19, 5, 12], the most extensive of which are those concerned with cardiac [20, 6, 28] and neuronal [8, 30, 33] cell dynamics. The corresponding equations for n electrically coupled FHN neurons are dxi dt = xi − xi 3 − yi + g n ∑ j=1 (xi − xj), dyi dt = (xi + a− b yi) (1) for i = 1, . . . , n and where we have taken a ∈ [−2, 2], b ∈ (0,∞), −1 ≤ g ≤ 1 (g 6= 0) and 0 < ≤ 0.1. It is well-known that in the FHN model, the variables have no direct physiological interpretation. However, for the parameter ranges quoted above, the qualitative behaviour of the x’s and y’s are similar to that of the voltage and gating variables in the HodgkinHuxley equations. †Correspondence to: IRCCyN, 1, rue de la Noe, BP 92101, 44321 Nantes Cedex 03 France. Tel : +33 2 40 37 69 00. Fax : +33 2 40 37 69 30 ar X iv :1 00 1. 54 20 v1 [ qbi o. N C ] 2 9 Ja n 20 10 Stationary points in the FHN system Neurons within the central nervous system can exist in a variety of dynamical states. Many, for example, are in a state of quiescence and elicit a relaxation oscillation, known as an action potential, when perturbed with a suprathreshold stimulus. This is oft termed as “spiking”. With the appropriate choice of parameters, the FHN model can exhibit such behaviour. However, upon varying these parameters, one can also invoke a Hopf bifurcation resulting in relaxation oscillations at an intrinsic frequency without the need for any stimulation. Such self sustained neurons can be found within the central nervous system and can have complex interactions with the environment. A typical example is that of circadian cells which act like the organisms clock cells and have been shown to have the ability to entrain their oscillations with the environments light-dark cycle [31]. Neurons have also been shown to express bistability [29],which again is present in the FHN model due to the cubic nonlinearity present in the system. The dynamical states in the FHN model described above can be identified via a stability analysis of the stationary points in the system corresponding to dxi dt = dyi dt = 0. The solutions to this system shall henceforth be referred to as simply the steady states. Such an analysis is thus essential if one is to have a basic understanding of the system. In fact, the steady state equations of this model have already been solved exactly for the n = 1 and n = 2 cases [9]. Recently a conjecture relating the steady states of a system and its Kolmogorov-Sinai entropy has brought to greater light the importance of steady state analysis for neural modeling in general[4, 1, 2, 3]. It should be noted that computational and numerical algebraic geometry methods have found many applications in many branches of theoretical physics in general (see, e.g., Refs. [24, 23, 14]). Below we use the DV method to solve the corresponding equations for the n = 3 case and study the stability and bifurcation structure of these systems. The functions we use are in the packages Groebner and RootFinding[Parametric] of the computer algebra system Maple 13. Due to the polynomial-like nonlinearity also exhibited in the HR, ML and Izhikevich models, the same technique may also be applied to these systems. The problem addressed in this paper can also be related to the more general problem of the algebraic analysis of the solutions of a differential system, studied for example in [7], [26],[32] or [34],[35]. It is also worth mentioning that, though we have performed our computations in Maple, some of the computations (quantifier elimination, sample points extraction) can be run using Discoverer∗, QEPCAD†, REDUCE‡, Mathematica§, STRINGVACUA¶ or RAG‖. Algebraic tools For this study, we used the Discriminant Variety [18] and Cylindrical Algebraic Decomposition [10]. The Discriminant Variety is an implicit representation of the desired partition, while Cylindrical Algebraic Decomposition describes explicitly each cell of the partition. Combined together, these two methods are well adapted to the analysis of the steady states (stable steady states). They provide a partition of the parameter space into con∗http://www.is.pku.edu.cn/ xbc/discoverer.html, developed by Wang, Xia et al. †http://www.usna.edu/Users/cs/qepcad/B/QEPCAD.html, developed by Hong et al. ‡http://reduce-algebra.sourceforge.net/, developed by Hearn, Weispfenning et al. §http://www.wolfram.com/products/mathematica/index.html ¶http://www-thphys.physics.ox.ac.uk/projects/Stringvacua/ ‖http://www-spiral.lip6.fr/ safey/RAGLib/distrib.html, developed by Safey El Din Figure 1: Parameter space decomposition adapted to the number of steady states for n = 3: in each connected component outside the surface, there is a constant number of steady states. g ] −∞, μ1[ ]μ1, μ2[ ]μ2, μ3[ ]μ3, μ4[ ]μ4, μ5[ ]μ5, μ6[ Steady states 3 3 3 3 3 3 Stable states 2 2 2 2 2 2 g ]μ6, μ7[ ]μ7, μ8[ ]μ8, μ9[ ]μ9, μ10[ ]μ10, μ11[]μ11, μ12[ Steady states 3 3 3 3 3 15 Stable states 2 2 2 2 2 2 g ]μ12, μ13[]μ13, μ14[]μ14, μ15[]μ15, μ16[]μ16, μ17[]μ17, μ18[ Steady states 15 15 15 15 27 27 Stable states 2 2 2 2 8 8 g ]μ18, μ19[]μ19, μ20[]μ20, μ21[]μ21, μ22[]μ22, μ23[]μ23, μ24[ Steady states 27 27 27 27 27 27 Stable states 8 8 8 14 6 6 g ]μ24, μ25[]μ25, μ26[]μ26, μ27[ ]μ27,∞[ Steady states 27 27 15 15 Stable states 6 6 6 6 where: μ1 ≈ −1.690122506, μ2 ≈ −.3010608116, μ3 ≈ −.2787081900, μ4 ≈ −.2530108085, μ5 ≈ −.2200916671, μ6 ≈ −.2119432516, μ7 ≈ −.1972651791, μ8 ≈ −.1944444444, μ9 ≈ −.1797585308, μ10 ≈ −.1768496368, μ11 ≈ −.1666666667, μ12 ≈ −.1111111111, μ13 ≈ −.09533095509, μ14 ≈ −.08765527212, μ15 ≈ −.08093319477, μ16 ≈ −.06753407911, μ17 ≈ −.06313592910, μ18 ≈ .1646239760, μ19 ≈ .2210124208, μ20 ≈ .2847525372, μ21 ≈ .3136620682, μ22 ≈ .3333333333, μ23 ≈ .3372096920, μ24 ≈ .3373316121, μ25 ≈ .3483417636, μ26 ≈ .3562092137, μ27 ≈ .3888888889. Figure 2: Number of steady and stable states for a = 0, b = 2 according to the parameter g. nected cells, such that within each cell, the number of steady states (stable steady states) is constant. Finally, the main algebraic criterion to decide if a solution is stable is the Routh-Hurwitz criterion, or its Liénard-Chipart variant [17]. For our problem, we used a reduced criterion, more adapted to the algebraic DV and CAD methods. Results The model in Eq. 1 has already been studied for the case n = 2 in Ref. [9]. We therefore focussed our work on the case n ≥ 3. First, we succeeded in describing the steady steady states for n = 3. The result of our computation is summarized in Figure 1. However, when n ≥ 3, the computations of the stable states are much more difficult and we did not succeed in describing the full parameter space according to the number of stable steady states. Nevertheless, by fixing the parameters a and b to the numerical values (respectively 0 and 2), we were able to describe the number of steady and stable steady states according to the free parameter g, for the n = 3 case. Since there is only one parameter, the description of the parameter space is a union of intervals in the parameter g for which the number of stable (and steady) solutions is constant. The results of the computations are summarized in Figure 2. Acknowlegement DM was supported by Science Foundation of Ireland. SP was in part supported by the Irish Research Council for Science Engineering and Technology (IRCSET). References and Notes [1] M. S. Baptista, F. M. Moukam Kakmeni, and C. Grebogi. The combined effect of chemical and electrical synapses in hindmarsch-rose neural networks on synchronisation and on the rate of information. arXiv:0910.0988v1 [nlin.CD], 2009. [2] M. S. Baptista, F. Moukam Kakmeni, Gianluigi Del Magno, and M. S. Hussein. The combined effect of chemical and electrical synapses in hindmarsch-rose neural networks on synchronisation and on the rate of information. arXiv:0805.3487v4 [nlin.CD], 2009. [3] M. S. Baptista and J. Kurths. Transmission of information in active networks. Phys. Rev. E, 77(2):026205, Feb 2008. [4] Murilo S. Baptista, Josue X. de Carvalho, and Mahir S. Hussein. Finding quasi-optimal network topologies for information transmission in active networks. PLoS ONE, 3:e3479, 2008. [5] S. Barland, O. Piro, M. Gludici, J.R. Tresicce, and S. Balle. Experimental evidence of van der polfitzhugh-nagumo dynamics in semiconductor optical amplifiers. Phys. Rev. E, 68(3):036209, 2003. [6] V. N. Biktashev, A. V. Holden, S. F. Mironov, A. M. Pertsov, and A. V. Zaitsev. Three-dimensional organisation of re-entrant propagation during experimental ventricular fibrillation. Chaos Sol. & Frac., 13(8):1713–1733, 2002. [7] F. Boulier, M. Lefranc, F. Lemaire, P.-E. Morant, and A. Ürgüplü. On proving the absence of oscillations in models of genetic circuits. In AB, pages 66–80, 2007. [8] David Brown, Jianfeng Geng, and Stuart Feerick. Variability of firing of hodgkin-huxley and fitzhughnagumo neurons with stochastic synaptic input. Phys. Rev. Lett., 82(23):4731–4734, 1999. [9] Sue Ann Campbell and Michael Waite. Multistability in coupled fitzhugh–nagumo oscillators. Nonlin. Anal., 47:1093–1104, 2001. [10] George E. Collins. Quantifier Elimination for Real Closed Fields by Cylindrical Algebraic Decomposition. Springer Verlag, 1975. [11] B. Van der Pol. A theory of the amplitude of free and forced triode vibrations. Radio Review, 1:701–710, 754–762, 1920. [12] Orazio Descalzi, Jaime Cisternas, Daniel Escaff, and Helmut R. Brand. Noise induces partial annihilation of colliding dissipitive solitons. Phys. Rev. Lett., 102(18):188304, 2009. [13] R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical J., 1(6):445–466, 1961. [14] James Gray, Yang-Hui He, Anton Ilderton, and Andre Lukas. STRINGVACUA: A Mathematica Package for Studying Vacuum Configurations in String Phenomenology. Comput. Phys. Commun., 180:107–119, 2009. [15] William Hanan, Dhagash Mehta, Sepanda Pouryahya, and Julien Sprott. Work in progress. [16] A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. (Lond.), 117:500–544, 1952. [17] P. Lancaster and M. Tismenetsky. The Theory of Matrices. Academic Press, New York, 1985. [18] Daniel Lazard and Fabrice Rouillier. Solving parametric polynomial systems. J. Symb. Comput., 42(6):636–667, 2007. [19] Anita Shu-Han Lin, Martin L. Buist, Nicolas P. Smith, and Andrew J. Pullan. Modelling slow wave activity in the small intestine. J. Theo. Bio., 242(2):356–362, 2006. [20] A.F.M Marée and A.V. Panfilov. Spiral breakup in excitable tissue due to lateral instability. Phys. Rev. Lett., 78(9):1819–1822, 1997. [21] F. Marino, M. De Rosa, and F. Marin. Canard orbits in fabry-perot cavities induced by radiation pressure and photothermal effects. Phys. Rev. E, 73:026217, 2006. [22] Bradley Marts, David J. W. Simpson, Aric Hagberg, and Anna L. Lin. Period doubling in a periodically forced belousov-zhabotinsky reaction. Phys. Rev. E, 76:026213, 2007. [23] Dhagash Mehta. Lattice vs. Continuum: Landau Gauge Fixing and ’t Hooft-Polyakov Monopoles. Ph.D. Thesis, 2009. [24] Dhagash Mehta, Andre Sternbeck, Lorenz von Smekal, and Anthony G Williams. Lattice Landau Gauge and Algebraic Geometry. POSQCD-TNT09:025, 2009. [25] J. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simulating nerve axon. Proc. IRE, 50(10):2061–2070, 1962. [26] Wei Niu and Dongming Wang. Algebraic approaches to stability analysis of biological systems. Mathematics in Computer Science, 1(3):507–539, 2008. [27] Hiroto Shoji, Kohtaro Yamada, and Takao Ohta. Interconnected turing patterns in three dimensions. Phys. Rev. E, 72:065202(R), 2005. [28] S. Sinha, K. M. Stein, and D. J. Christini. Critical role of inhomogeneities in pacing termination of cardiac reentry. Chaos, 12(3):893–902, 2002. [29] I. Tasaki. Demonstration of two stable states of the nerve membrane in potassium-rich media. J. Physiol. (Lond.), 148(2):306–331, 1959. [30] Rail Toral, C. Masoller, Claudio R. Mirasso, M. Ciszak, and O. Calvo. Characterization of the anticipated synchronization regime in the coupled fitzhugh-nagumo model for neurons. Physica A, 325(12):192–198, 2003. [31] Shobi Veleri, Christian Brandes, Charlotte Helfrich-F́’orster, Jeffrey C. Hall, and Ralf Stanewsky. A self-sustaining, light-entrainable circadian oscillator in the drosophila brain. Curr. Biol., 13(20):1758– 1767, 2003. [32] Dongming Wang and Bican Xia. Stability analysis of biological systems with real solution classification. In ISSAC ’05: Proceedings of the 2005 international symposium on Symbolic and algebraic computation, pages 354–361, New York, NY, USA, 2005. ACM. [33] D.Q. Wei, X.S. Luo, and Y.L. Zou. Firing activity of complex space-clampled fitzhugh-nagumo neural netowrks. Eur. Phys. J. B, 63:279–282, 2008. [34] H. Yoshida, K. Nakagawa, H. Anai, and K. Horimoto. An algebraic-numeric algorithm for the model selection in kinetic networks. In CASC, pages 433–447, 2007. [35] H. Yoshida, K. Nakagawa, H. Anai, and K. Horimoto. Exact parameter determination for parkinson’s disease diagnosis with pet using an algebraic approach. In AB, pages 110–124, 2007.